Sensitivity of a Shallow- Water Model to Parameters. 



Eugene Kazantsev 

INRIA, projet MOISE, Laboratoire Jean Kuntzmann, 
BP 53, 38041 Grenoble Cedex 9, France 
Telephone: +33 4 76 51 42 65 
Fax: +33 4 76 63 12 63 

E-Mail: kazan@iniag.fr 



Abstract: An adjoint based technique is applied to a shallow water model in order to 
estimate the influence of the model's parameters on the solution. Among parameters 
the bottom topography, initial conditions, boundary conditions on rigid boundaries, 
viscosity coefBcients Coriolis parameter and the amplitude of the wind stress tension 
are considered. Their influence is analyzed from three points of view: 

• flexibility of the model with respect to a parameter that is related to the low- 
est value of the cost function that can be obtained in the data assimilation 
experiment that controls this parameter; 

• possibility to improve the model by the parameter's control, i.e. whether the 
solution with the optimal parameter remains close to observations after the end 
of control; 

• sensitivity of the model solution to the parameter in a classical sense. That 
implies the analysis of the sensitivity estimates and their comparison with each 
other and with the local Lyapunov exponents that characterize the sensitivity 
of the model to initial conditions. 

Two configurations have been analyzed: an academic case of the model in a square box 
and a more realistic case simulating Black Sea currents. It is shown in both experiments 
that the boundary conditions near a rigid boundary influence the most the solution. 
This fact points out the necessity to identify optimal boundary approximation during 
a model development. 

Keywords: Variational Data Assimilation; Sensitivity to parameters; 
Boundary conditions; Shallow water model. 

1 Introduction 

Lorenz, in his pioneering work [I has shown that a geophysical fluid is ex- 
tremely sensitive to initial conditions and perturbations of the initial state grow 
exponentially in time. This discovery led to the development of the sensitivity 
studies intended to describe the evolution of a small unknown error in initial 
data and its influence on the forecast. 

Together with the sensitivity studies, data assimilation methods have also 
been under rapid development. These methods are intended to bring the model 
and various observational information together in order to better identify the 
initial state of the model. 
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The variational data assimilation technique, first proposed in [2] , [3] , is based 
on the optimal control methods [4] and perturbations theory [5] . This technique 
allows us to retrieve an optimal data for a given model from heterogeneous 
observation fields. Since the early 1990s, many mathematical and geophysical 
teams are involved in the development of the data assimilation strategy. One 
can cite many papers devoted to this problem, as in the domain of development 
of different methods for the data assimilation and also in the domain of its 
applications to the atmosphere and oceans. 

In the beginning, data assimilation methods were intended to identify and 
reconstruct an optimal initial state for the model. However, the idea that other 
model's parameters should also be identified by data assimilation has also been 
studied and discussed in numerous papers. One can cite several examples of 
using data assimilation to identify the bottom topography of simple models ([6], 
H, [H]); to control open boundary conditions in coastal and regional models ([S], 
Cn], [n], [H]), boundary conditions on rigid boundaries ([13], [H] [TS], [Hj, [H], 
[15] ) and to determine other parameters of a model ([H], [10], [21] )■ 

This paper is devoted to the comparison and the analysis of the dependence 
of the solution of a shallow-water model on different parameters. All available 
model's parameters have been included in the test set. First of all, among these 
parameters we consider the initial state of the model in order to compare the 
influence of all other model's parameters with the influence of the initial state. 
Second, we study the influence of the boundary conditions on rigid boundaries. 
But, instead of considering the boundary conditions themselves, we remind that 
particular attention must be paid to the discretization process (see [M]) that 
introduces the boundary conditions into the model. So, we use the representa- 
tion proposed in [16] that allows us to control the discretization of the model's 
operators near the boundary and to obtain immediately the numerical scheme 
that deal with the boundary conditions. Third, we also include the bottom 
topography in the set of control parameters because the importance of optimal 
representation of the topography on the model grid has been discussed in nu- 
merous papers, as cited above and others. In addition, the control set includes 
also several scalar parameters like dissipation coefficients and forcing amplitude. 
All these parameters have empirical origins, their values can also be optimized. 
And finally, we control the Coriolis parameter. 

The influence of the parameters on the model solution is analyzed from 
three points of view. First, we compare the flexibility of the model with respect 
to a particular parameter. The flexibility in this paper is understood as the 
capability of a parameter to bring the model's solution closer to observations 
reducing the cost function value in the assimilation procedure. Of course, obser- 
vational data contain measurement errors and result from the different physical 
processes in the nature. Consequently, there is no hope to get vanishing cost 
function varying any model's parameter. However, if variation of some param- 
eter allows to bring the model's solution closer to observations, we can consider 
the model is more flexible with respect to this parameter. 

However, flexibility of the model cannot be considered as a real improvement 
of the model's solution. Low cost function's value is obtained in frames of the 
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minimization procedure, under strong external control that force the model 
solution toward observational data. The distance between the solution and 
observations can rapidly increase after the end of the assimilation when the 
model is no longer forced to remain close to observations. That is why we 
analyze also the behavior of the solution beyond the assimilation window and 
examine the distance on longer time scales considering that data assimilation 
improves the model if the solution remains close to observations after the end 
of control. 

And the third sensitivity estimate considered in this paper is similar to the 

classical sensitivity characteristic that relates the norm of the perturbation of 
the solution and the norm of the perturbation of the parameter. Local (or 
finite time) Lyapunov exponents, for example, are evaluated using this ratio 
with initial conditions used as a parameter. This ratio shows how much a small 
normalized error in a parameter perturbs the solution at a given time. 

Two examples are considered in this paper: an academic case of the model in 
a square box with artificially generated observations and a more realistic case of 
assimilation of real observational data in the Black sea model. We should note 
that the same model is used in both configurations, but the difference consists 
in chaotic behavior of the model in the square box and regular flow in the Black 
sea. Temporal variability of the solution in the last case is only due to variations 
of the wind stress on the sea surface. 



2 Shallow Water Model 

2.1 Model's equations and discretization 

In this paper we consider a shallow-water model written in a conservative form: 



huv — fJ'h— 1 — ahu + tqTx 



dy\ dy 

^ = 

- —lhv'^+gh{h-H)-iJ.h—\-ahv + ToTy, (1) 

dh dhu dhv 

dt dx dy 

where hu{x,y,t) and hv{x,y,t) are two flux components that represent the 
product of the velocity by the ocean depth, h{x,y,t), that corresponds to the 
distance from the sea surface to the bottom of the ocean. The sea surface 
elevation is represented by the difference h{x,y,t) — H{x,y), where H{x,y) is 
the bottom topography. The model is driven by the surface wind stress with 
components Tx{x,y,t) and Ty{x,y,t) normalized by tq and subjected to the 
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bottom drag that is parametrized by linear terms ahu and ahv. Horizontal 
eddy diffusion is represented by harmonic operators div{^h'S/u) and div[^hS/v). 
Coriolis parameter is represented by the variable f{y) that is equal to /o + /3y 
assuming /3-plane approximation. Parameter g is the reduced gravity. The 
system is defined in some domain with characteristic size L requiring that 
both hu and hv vanish on the whole boundary of J7. No boundary conditions 
are prescribed for h. Initial conditions are defined for all variables: hu, hv and 
h. 

As usual, initial conditions are considered as the control parameter of the 
model in this paper. We study the sensitivity of the model to its initial point 
and assimilate data to find its optimal value. However, in addition to initial 
conditions, all other parameters of the model, and namely its bottom topography 
H{x,y), Coriolis parameter / = /(y), scalar coefficients ^,a,g and tq, are also 
considered as control variables. All of them are allowed to vary in the data 
assimilation procedure in order to bring them to their optimal values. The 
sensitivity of the model with respect to small variations of these parameters will 
also be studied. 

As it has been shown in [18], [17], the influence of the boundary conditions 
on rigid boundaries on the model's solution is also strong. On the one hand, 
assimilating data allows us to better represent the model's boundary on the 
model's grid, and, on the other hand, it allows to correct several types of errors 
committed by numerical schemes and finite dimensional approximations of the 
model's operators. 

However, following [17 , instead of controlling boundary conditions, we choose 
to control the way they are introduced into the model. And namely, we consider 
the discretization of the model's operators near the boundary as the control pa- 
rameter, that means the numerical scheme that takes into account the set of 
boundary conditions. A detailed description of controlling the discretization as 
well as the analysis of the tangent and adjoint models can be found in |18| . 
Here, we just remind the principal points. 

We discretize all variables of this equation on the regular Arakawa's C-grid 
[22] with constant grid step Sx — in both x and y directions: 

huij_i/2{t) = hu{ih,jh-h/2,t)ioTi = 0,...NJ = 0,...,N + l 
hvi_i/2j{t) = hv{ih- h/2,jh,t) ioT i ^ 0,. . .N + ^0,. ..N 
hi^i/2,j-i/2it) h{ih- h/2Jh- h/2,t) for i = 0,.. .N + l,j ^ 0,. ..N + 1 

Discretizing system (|T]), we replace the derivatives by their finite difference 
representations and Dy and introduce two interpolations in x and y coordi- 
nates Sx and Sy. Interpolations are necessary on the staggered grid to calculate 
the variable's values in nodes where other variables are defined. The discretized 
system (|T]) writes 

^ - fSxSyhv + Dx(^^^^+gh{h~H{x,y))-^^hDx^^ + 
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/ (Syhu Sxhv) , r. , hu \ 
+ Dy\^ S S h tJ'[SxSyh)Dy-^j = -ahu + TQTx, 

+ fSySxhu + D^l^ — ^ g, ^ ^.{,SxSyh)Dx-^ 

+ Dy( ^^y^^^ +gh{h-Hix,y))-^ihDypj-^ =-ahv + ToTy,{2) 



dh 
'dt 



-Dxhu — Dyhv. 



Discretized operators Dx,Dy and Sx,Sy are defined in a classical way at all 
internal points of the domain. For example, the second order derivative and the 
interpolation operator of the variable hu defined at corresponding points write 

{Dxhu)i_i/2j-i/2 ^ for I = 2,...,iV- 1, 

iSxhu)i^i/2j^i/2 = — ^ for I = 2,...,iV- 1. (3) 

To calculate fourth order approximations of derivatives and interpolations we 
use the following formulas 

,„ , - hui_2,j-i/2 - 27hui_i_j_i/2 + 27/iUi.j_i/2 - hui+ij^i/2 

{Dxhu),^l/2,j-l/2 = , 

,„ , ^ -hu,_2.j-i/2 + Qhui_i^j_i/2 + 9huij_i/2-hui+ij_i/n 
[Sxnu)i^i/2.j-i/2 = 

Discretization of operators in the directly adjacent to the boundary nodes 
are different from and represent the control variables in this study. In order 
to obtain their optimal values assimilating external data, we suppose nothing 
about derivatives and interpolations near the boundary and write them in a 
general form 



Sx 



D 

{Dxhu)i/2j-i/2 = "o " 

[Sxhu)i/2J-l/2 ^ 2 

This formula represents a linear combination of values of hu at two points 
adjacent to the boundary with coefficients a. The constant ao may be added in 
some cases to simulate non-uniform boundary conditions like hu{0, y) = ^ 0. 

We distinguish a for different variables and different operators allowing dif- 
ferent controls of derivatives because of the different nature of these variables 
and different boundary conditions prescribed for them. It is obvious, for ex- 
ample, that the approximation of the derivative Dx in the first equation may 
differ from the approximation of in the third one. Although both oper- 
ators represent a derivative, boundary conditions for hu and h are different; 
these derivatives are defined at different points, at different distance from the 
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boundary. Consequently, it is reasonable to let them be controlled separately 
and to assume that their optimal approximation may be different with distinct 
coefficients a^== and a^^ . 

Time stepping of this model is performed by the leap-frog scheme. The first 
time step is split into two Runge-Kutta stages in order to ensure the second 
order approximation. 

The approximation of the derivative introduced by (O and (O depends on 
variables a. They are added to the set of control variables enumerated above. 
Operators are allowed to change their properties near boundaries in order to 
find the best fit with requirements of the model and data. To assign all control 
variables we shall perform the data assimilation procedure and find their optimal 
values. Variational data assimilation is usually performed by minimization of 
the specially introduced cost function. The minimization is achieved using the 
gradient of the cost function that is usually determined by the run of the adjoint 
to the tangent linear model. 

2.2 Cost function 

One of the principal purposes of variational data assimilation consists in the 

variation of control parameters in order to bring the model's solution closer 

to the observational data. This implies the necessity to measure the distance 

between the trajectory of the model and the data. Introducing the cost function, 

we define this measure. Generally speaking, the cost function is represented by 

some norm of the difference between the model's solution and observations, 

eventually accompanied by some regularization term. 

To define the cost function we introduce dimensionless state vector that is 

composed of three variables of the model — {whuhu, Wkvhv, WhhY weighted by 

coefficients w. These weights are used to normalize values of the flux components 

by Whu = Whv — T? — 1 and the Sea surface elevation by Wh = -n-- The 
Hoy/gHo ^0 

distance between the model solution and observations is defined as the Euclidean 
norm of the difference 



In this expression, we emphasize implicit dependence of ^ on time and on the 
set of the control parameters p that is composed of 

• the set of initial conditions of the model (pQ — {hu \t=o, hv \t=Q, h |t=o}j 

• the set of the coefficients a that controls the discretizations of operators 
near the boundary, 

• the bottom topography H{x, y) 




(6) 



k 
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• four scalar parameters a, fi, g, tq 



• and the Coriolis parameter /(y). 
Taking into account the results obtained in ^17j . we define the cost function as 

T 

Ap) - / te{HP,t))dt (8) 



that gives bigger importance to the difference f ^ at the end of assimilation 
interval. 

It should be noted here, that this cost function can only be used in the 
case of assimilation of a perfect artificially generated data. When we assimilate 
some kind of real data that contains errors of measurements and is defined on 
a different grid, we should add some regularization term to the cost function 
(like the distance from the initial guess) and use some more appropriate norm 
instead of the Euclidean one (see, for example [53] for details) . 

The nth component of the gradient of the cost function can be calculated as 
the Gateaux derivative of an implicit function: 



T T 

t2\ r / «C2 



Ok 





T 



dt 



dpn 

2/t(E(^.-0f^)|^)rf^ (9) 



because the derivative can easily be calculated from — 2( 

'). The second term in (l9|), jv^, represents the matrix of the tangent linear 

model that relates the perturbation of the parameter p„ and the perturbation 
of fcth component of the model state vector (t>k. This relationship, of course, 
is assumed in the linear approach, that means it is only valid for infinitesimal 
perturbations. 

In the classical case, when initial conditions are considered as the only con- 
trol variable, the derivative ^'q^^ = '^dcj)^ classical tangent model that 
describes the temporal evolution of a small error in the initial model state. 
The matrix is a square matrix that is widely studied in numerous sensitivity 
analyses. Its singular values at infinite time limit are related to well known 
Lyapunov exponents that determine the model behavior (chaotic or regular) 
and the dimension of it's attractor. 

In our case, the matrix '^^^ is rectangular. It describes the evolution of an 

infinitesimal error in any parameter (including initial state). However, we can 
study its properties in a similar way as we do with the classical tangent linear 
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model. Its structure and composition is described in [T7] for the case of using 
coefficients a as control parameters and in [8] for the case when the bottom 
topography is used to control the model's solution. 

The product X]fc(<^fc ~ '^fc'"')^^ ® represents an unusual vector- matrix 
product. To calculate this product directly we would have to evaluate all the 
elements of the matrix. This would require as many tangent model runs as the 
size of the state vector is. So, instead of the tangent model, we shall use the 
adjoint one that allows us to get the result by one run of the model. Backward in 
time adjoint model integration that starts from ((/) — 0°''*) provides immediately 



the product " <^°'") which is exactly equal to - (f^')^ in ©. 

Using these notations, we write 

T ^ 

VI = 2 1 i (^^^ i<t>ip, t) - r'' {t))dt (10) 



where the expression in the integral is the result of the adjoint model run from 
t to starting from the vector {(j){p,t) — 

Tangent and adjoint models have been automatically generated by the Tape- 
nade software [25, [25] developed by the TROPICS team in INRIA. This soft- 
ware analyzes the source code of the nonlinear model and produces codes of its 



derivative -rr- and of the adjoint 

op \.op 

This gradient is used in the minimization procedure that is implemented in 

order to find the minimum of the cost function: 



I{p) = minl(p) (11) 

p 

Coefficients p are considered as coefficients achieving an optimal parameter for 
the model. As it has been already noted, the set of parameters p is composed of 
the set of initial conditions of the model (j)o , the set of the coefficients a that con- 
trols the discretization of operators near the boundary, the bottom topography 
H{x,y) four scalar parameters a,fj,,g,To and the Coriolis parameter f{y). We 
shall minimize the cost function controlling as the total set of available param- 
eters p and any possible subset comparing the efficiency of the minimization. 

We use the minimization procedure developed by Jean Charles Gilbert and 
Claude Lemarechal, INRIA ^ . The procedure uses the limited memory quasi- 
Newton method. 



2.3 Sensitivity estimates 

In addition to the data assimilation, we perform also the sensitivity study of the 
model solution to parameters enumerated in the previous subsection. We are 
looking for a perturbation in the model's parameters 6p that, for a given small 
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norm, maximizes the norm of the perturbation of the sohition at time t. 

We cannote that we aheady have all the necessary software to estimate X{t). 
Tangent linear model allows us to calculate S(f>{t) — (^~^^^ ^P- Us- 

ing the scalar product that corresponds to the norm in the definition of the 
distance ^ (0), we can write 

mt)\. (dm 



\{t) = max ' — max 



<.dp, 5p:> <:5p, 6p:> 



= max ; — ; (13) 

This expression is a well known Rayleigh-Ritz ratio which is equal to the largest 
eigenvalue of the problem 

So far, we need just the maximal eigenvalue and the matrix of the problem is a 
self-adjoint positive definite matrix, we can solve the problem (|14l) by the power 
method performing successive iterations 



(dmy (debit) 

I dp J I dp 



dn+i = — 7 ( * ) T , ^0 = random vector 

In the limit, the denominator of the right- hand-side tends to the largest eigen- 
value and to the corresponding eigenvector of the matrix. The principal 
advantage of this method consists in the fact that we do not need to calculate 
the matrix itself, we just need a matrix-vector product. So far, we have both 
codes for the tangent and adjoint models; we can successively run these models 
and get the left-hand side of (|14|) . 

We should note here that when the initial conditions of the model are used as 
the control parameters (i.e. 5p = d4>{0)), the sensitivity characteristics X{t) are 
all close to one when t — >■ 0. It is evident because the perturbation has no time 
to be transformed by the model's dynamics and we get S(j){t) \t ^q— (50(0) = dp. 

When any other model parameter is used as the control and the error growing 
time is small, all A(t) are vanishing. This is also clear: the model's dynamics 
has no time to transmit the perturbations from the parameters to the solution. 
The perturbation of the solution remains, consequently, close to zero as well as 
the value of X{t) \t >o= 0. 



9 



In order to make the behavior of the sensitivity characteristics uniform with 
different parameters, we shall use X{t) — 1 every time when the initial model's 
state is considered as the control parameter. 

Another point that we should emphasize, concerns the physical dimensions 
of parameters. So far we want to compare the sensibility of the model to vari- 
ous parameters; we should be careful with bringing them to the dimensionless 
form because the choice of the characteristic values influences the result. In 
this paper, we choose to measure all the perturbations in fractions of the orig- 
inal non-perturbed parameter. This will ensure relatively uniform weighting 
of perturbations. That means the perturbation of the Coriolis parameter is 
normalized by /o, the value of this parameter in the middle of the basin; the 
perturbation of the bottom topography is normalized by Hq, the average depth; 
perturbations Sfi, Sa, 5g and Stq are respectively normalized by /x, a, g and 
tq. As it has been already noted, perturbations of the initial conditions i5^o 
and of the model state vector 5(t){t) are already dimensionless, being obtained 

form the model's variables using weights Whu = w^v — tt — \ tt and Wh = -n-- 

HoVgHo ^0 

Coefficients a that are used to control the discretization of operators are also 
already dimensionless, having characteristic values around one (-1-1 or -1 in the 
second order derivatives and 1/2 in the interpolations) they are used without 
normalization. 

Thus, the Rayleigh-Ritz ratio ([T^ and the eigenvalues of problem ([Tl)) be- 
come dimensionless, but they keep the dependence on the normalization con- 
stants. 



3 Model in a square box. 

We start from the data assimilation in frames of the very well studied "aca- 
demic" configuration. Several experiments have been performed with the model 
in a square box of side length L — 2000 km driven by a steady, zonal wind 
forcing with a classical sinusoidal profile 

27r(u-i/2) 

Tx = Tq cos 

that leads to the formation of a double gyre circulation [57] . The attractor of the 

model and the bifurcation diagram in a similar configuration has been described 

in [28]. Following their results, we intentionally chose the model's parameters to 

ensure chaotic behavior. The maximal wind tension on the surface is taken to 

be To = 0.5 , The coefficient of Eckman dissipation and the lateral friction 
cm 

2 

coefficient are chosen as cr = 5 x lO^^s^^ and fi — 200^^ respectively. 

As it has been already noted, the Coriolis parameter is a linear function in 
y with fo = 7x 10~^s~^ and (3 = 2 x 10~"'^^(tos)"^. The reduced gravity and 
the depth are respectively equal to g — 0.02^, Hq — lOOOm. 

The resolution of the model in this section is intentionally chosen to be too 
coarse to resolve the Munk layer [29] that is characterized by the local equilib- 
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rium between the /3-efFect and the lateral dissipation. Its characteristic width is 



present case. The model's grid is composed of 30 nodes in each direction, that 
means the grid-step is equal to 67 km, that is more than the Munk parameter. 
Thus, there is only one grid node in the layer and the solution exhibits spurious 
oscillations near the western boundary due to unresolved boundary layer. 

Artificial "observational" data are generated by the same model with all the 
same parameters but with 9 times finer resolution (7.6 km grid step). The fine 
resolution model, having 7 nodes in the Munk layer, resolves explicitly the layer 
and must have no spurious oscillations. All nodes of the coarse grid belong to 
the fine grid, consequently, we do not need to interpolate "observational" data 
to the coarse grid. We just take values in nodes of the high resolution grid that 
correspond to nodes on the coarse grid. 

The model on the fine grid has been spun up from the rest state during 3 
years. The end of spin up was used as the initial state for the further integration 
of the model. From the result of this integration we have extracted values of 
all three variables at all grid points that belong to the coarse grid (as it has 
been noted, the grids have been chosen so, that all grid points of the coarse 
grid belong to the fine grid). This set is used as artificial observations in the 
following experiments. 

So far the model is nonlinear with intrinsically instable solution, there is no 
hope to obtain close solutions in long time model runs because any difference 
(even infinitesimal) between two models grows exponentially in time. Conse- 
quently, we have to confine our study to the analysis of a short time evolution 
of the model's solution simulating the forecasting properties of the model. 

All operators in the model are approximated either with the second or with 
the fourth order accuracy in the interior of the domain by ([3]) or (|4]). Initial 
guess for the coefficients a is defined to satisfy the second order scheme. That 
means the expression ([5]), that is used to interpolate functions and to calculate 
their derivatives near the boundary, is written with ao = 0. Coefficients a are 
defined as af = — 1, = 1 for all derivative operators and af = af = 1/2 
for all interpolations. That gives, for example, the value of the derivative of u 

at the point i = 1/2 as (-D£cm)i/2j"-i/2 = ^'^^^^"^ 6x 
3.1 Data assimilation 

Final values of the cost function obtained in experiments of the assimilation 
of artificially generated data into the shallow water model are shown in the 
table 1. All the experiments are carried out in the same conditions with the 
same data. The assimilation window has been chosen as 5 days interval. This 
time interval corresponds well to the characteristic time of the model's physics. 
Gravity waves, with their velocity equal to y/gHo = 4.5^, cross the 2000 km 
box in 5 days. 

As the initial guess for the initial conditions we use the state vector of the 




1/3 



which is equal to 42 km in the 
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high resolution model reduced on the coarse grid. This state is also used as the 
initial conditions in all other assimilation experiments with other control param- 
eters. Noted above values of the model's parameters (flat bottom topography, 
linear in y Coriolis parameter and scalar parameters (p, a, tq, g) arc used as 
the initial guess in the experiments that control these parameter; otherwise, we 
simply use these parameters in the model. 

As it has been noted, we control the parameters in the second and the fourth 
order model. With no assimilation at all, the solution of the fourth order model 
is "worse" than the solution of the second order one. That means the model's 
trajectory moves away from observations more rapidly and produces bigger cost 
function value. The reason of this is clear: high order scheme works worse 
when principal physical scales are not resolved explicitly. Indeed, the coarse 
grid of the model does not resolve the Munk boundary layer. The grid step 
5x = 67 km is bigger than the Munk parameter d = 2(/i//?)^/'^ = 42 km. But, it 



is the ratio ^ j where n is the order of approximation that determines the 



approximation error. Higher the order of approximation, bigger is this ratio and 
bigger is the error in the approximation of the boundary layer. Consequently, it 
is more important to identify an optimal numerical scheme for approximation 
of the boundary layer for the fourth order model. Indeed, if we assimilate data 
and find optimal parameters of the model, the fourth order model becomes 
comparable and even "better" (allowing lower cost function's value) than the 
second order one. 

However, the influence of different parameters on the solution is not the 

same. Comparing the cost function's values at the end of minimization proce- 
dure, we can say that the model is the most flexible with respect to the control 
of initial and boundary conditions. The cost function value can be divided by 3 
or even 4 in the minimization procedure. Controlling all other parameters, we 
cannot achieve such a low cost function. We must note also that the control of 
coeflacients a is much more expensive than the control of initial conditions: 5 
and even 10 times more iterations are required for the minimization to converge. 

Obviously, the best result is obtained in the joint control of all available 
parameters of the model. The cost function value is divided by 10 but the 
number of iterations exceeds 100. 

Along with the value of the cost function in the assimilation window, we 
examine the behavior of the model beyond the window. In fact, assimilating 
known data we can find optimal model parameters, but the optimality is guar- 
anteed in the assimilation window only. However, if the model is developed to 
provide a forecast, it has to use the observational data in the past and predict 
the future. In the described test case, consequently, it is more interesting to 
see the behavior of the solution after the end of assimilation, i.e. beyond the 
assimilation window. 

In order to see whether optimal parameters can improve the model's behavior 
even when the assimilation is over, we examine the difference £^{t) between the 
solution and observations over 20 days interval, i.e. during 15 days next to the 
assimilation window. The evolution of the difference ^{t) for different optimal 
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Control 


Second order interior 
Ifinai Iterations 


Fourth order interior 
^finai Iterations 


Nothing 


2.51 




3.29 




Initial state (f>Q 


1.04 


14 


0.96 


13 


Topography H{x,y) 


1.81 


5 


2.37 


7 


Coriohs f{y) 


2.15 


57 


2.81 


54 


Scalars cr, n,To,g 


2.20 


33 


3.06 


31 


Boundary a 


1.41 


113 


0.95 


124 


All parameters together 


0.34 


115 


0.28 


149 



Table 1: Number of iterations and final values of the cost function in the data 
assimilation experiments with different control parameters 

subsets of the controllable set is shown in figHJ Upper solid line in this figure 
represents the distance "observations-model" with no assimilation at all. Default 
parameters, described above, are used in this model run. 



Distance Distance 




Figure 1: Evolution of the distance £,{t) ^ during and after assimilation 
obtained with the second order model (left) and the fourth order model (right). 

Analysing the figure figlU we cannote several differences in the influence of 
model parameters on the solution. As it has been already seen, the flexibility 
of the model is small with respect to the bottom topography, scalar coefficients 
and the Coriolis parameter (joint control of these two parameters is presented). 
As a consequence, in the assimilation window (0 < t < 5 days) corresponding 
lines are close to the original line obtained with no assimilation at all. Beyond 
the window, these lines remain parallel to the original line keeping the difference 
that has been obtained in the minimization. In some cases, the increase of the 
distance from observations can even be less rapid than the increase of the solid 
line. This fact indicates that the optimal values of parameters obtained in the 
assimilation window remain to be optimal even beyond the window allowing us 
to use them for the forecast. 
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On the other hand, the model started from the optimal initial state (po 
exhibits a very different behavior in the window ensuring low cost function value. 
However, just after assimilation end, the distance from observations starts to 
increase and moves toward the solid line obtained with the default parameters. 
That means, controlling the initial state of the model allows us to improve the 
model in the assimilation window and the short range forecast (5-10 days in this 
experiment) as well, but has almost no influence on longer forecasts (15 days 
here). 

The most spectacular result of the model improvement can be obtained in 
the experiment that controls coefhcients a for the fourth order model. The value 
of the distance from observations is divided by 2 at T = 20 days. That means 
optimal parametrization of boundary conditions remains optimal after the end 
of assimilations and can help to improve the model bringing the solution closer 
to the solution of a high resolution model used to generate artificial observa- 
tions in this experiment. This can be explained by the improved accuracy of 
the fourth order approximation in the interior of the domain accompanied by 
optimal boundary conditions that are really necessary for this model. 

Thus, we see in this experiment that if the model's dynamics suffers from low 
resolution and other numerical errors, better forecast is achieved by controlling 
the model's operators rather than initial conditions. 

3.2 Sensitivity estimates. 

The third way of the sensitivity analysis consists in solving of the eigenvalue 
problem (|14p and analyzing X{t) on different scales of error growing time from 
about 10 minutes (10^^ day) to more than one year (500 days). Characteristic 
time scale (5 days during which the gravity waves cross the domain) is situated 
in the middle of this interval. Performed experiments with the second and the 
fourth order models show that X{t) are very close to each other. So, we plot 
only one of them in figlU and namely sensitivity estimates of the second order 
model. As it has been already noted, X{t) — 1 is plotted in the case when initial 
conditions are considered as the parameter. 

Analyzing the figure figlH we can see that three time scales can be clearly 
distinguished for the sensitivity characteristics of the model. The first, short 
time scales, approximately from up to 2-3 hours is characterized by the lin- 
ear growth of X{t). We should note that the growth is not only linear in the 
logarithmic coordinates, but the slope is equal to 1. Consequently, X(t) is pro- 
portional to t. Indeed, the model behaves as a linear model on these scales, the 
model's solution can be well approximated by just one step of the numerical 
time scheme. 

The second time scale that can be distinguished in the figure figH] corre- 
sponds to error growing times from 2-3 hours to 10-100 days. On these time 
scales we see slower growth of the sensitivity characteristics X{t) and, sometimes, 
no growth at all. These time scales are characterized by the modification of the 
stable-instable subspaces of the model. Instable space on short time scale is not 
the same as for long time evolution. Short time instabilities are usually localized 
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Singular values Singular values 




Figure 2: Sensitivity characteristics X{t) as functions of the error growing time 
(Log-Log coordinates (left) and Log-Linear coordinates (right). 

in space, while long time eigenvectors of (fT4|) possesses a global structure. On 
the medium time scales we see the transformation of the instable space of the 
quasi-linear model to the instable space of the non-linear chaotic model. This 
modification stipulates the slowdown (and even stagnation) of growth of A(t). 

The third time scale corresponds to the error growing times more than 100 
days. On these scales the model exhibits non-linear chaotic behavior with ex- 
ponential growth of all X{t). In order to zoom these time scales, we plot the 
same data in the Log-Linear coordinates in figI2]on the right. One can see that 
the growth on this time scale is purely exponential with the same exponent 
X{t) — ^exp(0.027t). The multiplier A is particular for each parameter, but 
the exponent is always the same. This confirms the remark made in [8], |17] : 
no matter how the perturbation was introduced into the model, its long-time 
growth is determined by the model's dynamics. 

Comparing the evolution of the sensitivity of the model to different param- 
eters, we see that on small scales it is the bottom topography that the model 
is the most sensitive to (thin solid line in figl2]). An error in the topography 
produces 13 times bigger perturbation in the model state than a similar error 
in the model's initial conditions (thick solid line in figEJ. However, X{t) does 
not grow at all on medium scales due to significant changes in the eigenvector's 
pattern. This leads to the fact that on long scales, the sensitivity of the model 
to the bottom topography is about 2 times lower than the sensitivity to initial 
conditions. 

On the other hand, the sensitivity of the model to the discretization of 
operators near the boundary exhibits the opposite behavior. On short scales, 
corresponding A is 2 times lower than A obtained for perturbations of 00, but 
there is no stagnation of the growth on the middle scales. As a result, we see 
that the model is 4 times more sensitive to a than to for long error growing 
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times. 



4 Model of the Black Sea. 

In this section we use the same model, but all the parameters are defined to 
describe the upper layer circulation of the Black sea. Configuration of the model 
and observational data have been kindly provided by Gennady Korotaev from 
the Marine Hydrophysical Institute, National Academy of Sciences of Ukraine, 
Sevastopol, Ukraine. This configuration is described in [50] . 

The model grid counts 141 x 88 nodes that corresponds to the grid box of 
dimension 7860 m and 6950 m in a; and y directions, respectively. A 15 minutes 
time step is used for integration of the model. The Coriolis parameter is equal 
to /o = lO^^'s^^ and /? = 2 x lO^^^(ms)^^. Horizontal viscosity is taken as 
= 50m'^s^^. Using a typical density difference between upper and underlying 
layers of 3.1kg/m^ , and unperturbed layer thickness of Hq = 150m, the Rossby 
radius of deformation is estimated at about 22 km and the reduced gravity value 
g = 0.031m/s^. The grid therefore resolves the mesoscale processes reasonably 
well. 

The model has been forced by the ECMWF wind stress data, available as 
daily averages for the years 1988 through 1999. Dynamical sea level recon- 
structed in |5T] was used as observational data in this section. These data have 
been collected in ERS-1 and TOPEX/Poseidon missions and preprocessed by 
the NASA Ocean Altimeter Pathfinder Project, Goddard Space Flight Center. 
Observational data are available from the 1st May 1992 until 1999. These data 
have been linearly interpolated to the model grid. 

So far the sea surface elevation is the only observational variable available 
in this experiment, we put Whu = Whv = in ([7]). Consequently, the difference 
between the model's solution and observations is calculated taking into account 
the variable h only. 

As it has been already noted, absence of observational data for the velocity 
fields brings us to modify the cost function. We have to add the background term 
in the cost function in order to require the velocity field to be sufficiently smooth. 
Otherwise, lack of information about velocity components in observational data 
would result in a spuriously irregular field obtained in assimilation. To ensure 
necessary regularity of hu and hv we add the distance from the initial guess 
to the cost function ([5]). In order to emphasize the requirement of smoothness, 
this distance is measured as an enstrophy of the difference between the initial 
guess and the current state: 



E 



dx dy 



(15) 



where hvP , hv^ denote flux components of the initial guess of the minimization 
procedure. 

Moreover, using real observational data requires to add at least one another 
term to the cost function. One can see that in the Figure 2 in [HP, spatially 
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averaged sea surface elevation of the Black sea exhibits a well distinguished 
seasonal cycle. That means the mass is not constant during a year, it decreases 
in autumn and increases in spring. Consequently, if we assimilate data during 
a short time (a season or less), we assimilate also the information about the 
mass flux specific for this season. This flux cannot be corrected later by the 
model because the discretization of operators near the boundary (that controls 
the mass evolution) is obtained once for all seasons. The mass variation of the 
Black sea reaches 25 centimeters of the sea surface elevation. Assimilating data 
within one season may, consequently, result in a persisting increase or decrease 
of the sea level of order of 50 cm per year. To avoid this spurious change of the 
total mass, we must either take the assimilation window of at least one year, or 
prescribe the mass conservation to the model's scheme. One year assimilation 
window is computationally expensive and is not justified by the model's physics. 
On the other hand, prescribed mass conservation removes just the sinusoidal 
seasonal variation, allowing us to keep all other processes and to choose any 
assimilation window we need. 

To correct the mass flux of the model, we add the following term to the cost 
function 




Similar to ([T5]), this term also ensures the regularity of the solution. It can 
be noted here that other terms may be added to the cost function in order to 
make a numerical scheme energy and/or enstrophy conserving, but we do not 
use them in this paper. 

The total cost function in this section is composed of three parts: ([8]), (|15l) 
and (US]): 

-^total — ^ ~\~ '^fl^srnooth + "f2^mass (17) 

Coefficients 7 are introduced to weight the information that comes from obser- 
vational data (with I) and an a priori knowledge about mass conservation and 
regularity of the solution. 

This modification of the cost function results, of course, in additional terms 
in the gradient: 

Vltotai = VI+271 foiDy {hu-huo)+DlD, {hv-hvoyj +272 J2 ("tyiu.j (0) 

(18) 

The model is spun up from the beginning of 1988 to May 1992 using the 
wind tension data on the surface. The state corresponding to the 1st of May 
1992 12h GMT is used as the initial guess in the data assimilation procedure 
controlling initial conditions of the model. The assimilation controls the initial 
conditions (pQ only with the assimilation window T — 1 day and the regular- 
ization parameter 71 = 0.04. Such a short window was chosen in order to get 
almost instantaneous state of the model to be used in further experiment as an 
initial state. 
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In this paper we chose a T = 30 days window which is longer than synoptic 
time scales. The minimization of the cost function has been accompanied by 
the mass preserving correction (jl6p with 72 = 0.01. 

4.1 Data Assimilation 

Like in the previous section, we examine the influence of the model parameters 
from two points of view. First, we assimilate observational data observing the 
cost function value in the assimilation window. This experiment reveals the 
flexibility of the model with respect to a parameter; it shows how close the 
solution can be with the optimal values of this parameter. Analysis of the cost 
function beyond the window shows the capability of the parameter to improve 
the forecast quality. And second, we apply the classical sensitivity analysis of 
the model solution with respect to parameters calculating the largest eigenvalue 
A(i) of problem (jl4[) that shows the influence of a small error in parameter on 
the solution of the model at time t. 

As we have already noted, we consider the same model, but there is a prin- 
cipal difference with the previous case of the model in a square box. The behav- 
ior of the model solution is not chaotic in this configuration. Variability of the 
model is generated directly by the variability of the wind stress on the surface. 
Consequently, we can compare particular trajectories of the model on any time 
interval because their evolution is stable without exponential divergence. Thus, 
we can hope that assimilating data in a relatively short window allows us to 
bring the model's solution closer to observation for a long integration period. 

The flexibility of the model is illustrated in figl3]on the left. We perform the 
data assimilation experiment with 30 days assimilation window using parameters 
described above as initial guess. Due to high CPU time of the data assimilation, 
we limit the number of iterations of the minimization procedure by 20. Thus, 
we have similar and reasonable computational cost in each experiment. 

One can see the model is the least flexible when we control the bottom 
topography and the most flexible with respect to the control of coefficients a. We 
can get 2 times lower distance solution-observations controlling discretization 
of operators near the boundary than controlling the topography. One cannote 
that Black sea model is less flexible with respect to the initial conditions than 
the model in the square box. This fact is due to the additional regularization 
term ([T5|) that is added to compensate the lack of observational data for u and v 
variables. On the other hand we can see increased flexibility with respect to the 
boundary conditions despite the presence of the forced mass conservation ()16p . 
This fact indicates the importance of controlling a for a model in a realistic 
domain. 

The evolution of the distance "model-observations" ^{t) is shown in the figl3] 
on the right. We assimilate the data in the 1 month window and integrate the 
model for the next 11 months. Assimilation window is zoomed in this figure in 
order to better distinguish different lines. One can see several differences with 
the behavior of the model in the square box. Controlling the initial state and 
the topography drastically increases the distance at the origin ^(0) indicating 
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Figure 3: Convergence of the cost function in the minimization procedure (left) 
and evolution of the distance "model-observations during 1 year (right). 



the initial guess is far from being optimal. The behavior of the model's solution 
with optimal initial point, with optimal topography and with optimal scalar 
parameters (not shown in the figure because the line is indistinguishable from 
other lines) are very similar as in the window and beyond. 

There are also common points with the experiment in the square box. It 
is the control of coefficients a that allows us to improve a long range forecast. 
Optimal initial point for the 30 days assimilation window influences the fore- 
cast during 100 days and after that there is no significant difference between 
the model with optimal parameters and the model with the default ones. An- 
other common point consists in an expected fact that the joint control of all 
parameters ensures the lowest distance from observations as in the window and 
beyond. 

4.2 Sensitivity estimates. 

And finally, we consider the sensitivity characteristics X{t) of the model to its 
parameters. The dependence of X{t) on the Error Growing Time (EGT) is 
shown in the figlH The figure is also divided into two parts: short time scales 
are zoomed on the left and long time scales on the right. Comparing this figure 
with the fig 12] we see that there is no general exponential growth of X{t) on long 
time scales. Moreover, the perturbation of initial conditions decreases on long 
time scales and corresponding X{t) become smaller than one (that is why we 
cannot plot X{t) — 1 in logarithmic coordinates for t > 20). This fact shows 
that the model solution is not chaotic and the variability of the solution is 
only due to the variability of the wind stress on the surface. The sensitivity 
of the model to other parameters does not show any regular behavior on long 
time scales. While A(t) that correspond to the boundary conditions and to 
the Coriolis parameter are growing with time (but not really exponentially), 
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the sensitivity to scalar model parameters oscillates and the sensitivity to the 
bottom topography stagnates. 

However, comparing fig|3] and figH] one can remark also several common 
points. We can also distinguish short, medium and long scales on which the 
behavior of X{t) is different. Obviously, linear error growth is also observed 
on short scales. On these scales, the model solution is also the most sensitive 
to topography perturbations and the sensitivity to boundary conditions is also 
smaller than the sensitivity to initial conditions. On the medium time scales we 
can also see common points with the square box. Sensitivity to the topography 
reaches the value A(i) — 2 and stagnates after that; A(t) corresponding to the 
boundary conditions and the Coriolis parameter continue to grow. 




Figure 4: Sensitivity characteristics X{t) as functions of the Error Growing Time 
(Log-Log coordinates (left) and Log-Linear coordinates (right). 



5 Conclusion 

The comparative study presented in this paper shows the influence of different 
model parameters on the solution. We do not show optimal parameters of 
the model that have been obtained in the assimilation process, nor the most 
sensitive patterns obtained as singular vectors of ([T4|) reserving all these data 
and analysis to a more technical study and discussion. In this paper we do just a 
comparison of the sensitivity of the model to the set of its internal and external 
parameters. 

The study is confined to the analysis of a low resolution model with a rather 
limited physics. Consequently we must acknowledge that the results may be 
valid only in the described case. Additional physical processes (baroclinic dy- 
namics, variable density due to heat and salinity fiuxes, etc.) may modify results 
of this study revealing other parameters the model may be sensitive to. 
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The main conclusion we can made from this comparison is the important role 
played by the boundary conditions on rigid boundaries. Ahuost all experiments 
show that the model is the most flexible with respect to control of coefficients 
a, this control allows us to bring the model's solution closer to the solution of 
the high-resolution model or to the observed data. In the experiment with the 
fourth order model in the square box, it is the control of the numerical scheme 
near the boundary that represents the major part of the total flexibility of the 
model (see fig [T] on the right). 

Optimal a found in the assimilation window remain optimal long time after 
the end of assimilation improving the forecasting ability of the model. We could 
see that the fourth order model in the square box allows us to divide by two 
the forecast error of the 20 days forecast. Optimal a obtained in one month 
assimilation remains optimal even for a one year run of the Black sea model. 

Finally, the long time sensitivity of the model's solution to a exceeds the 
sensitivity to almost all other parameters including the sensitivity to initial 
conditions. A perturbation of a of a given small norm results in a bigger per- 
turbation of the model's solution than a perturbation of some other parameter 
of an equal norm. 

However, we could see that the influence of boundary conditions is only 
important on long time scales, i.e. time scales that exceeds the characteris- 
tic time of the domain. In both experiments presented above the character- 
istic time was approximately equal to 5 days (as it was mentioned above, we 
use Tchar = L/^/gHo), and in both experiments the sensitivity to a becomes 
important on scale longer than 5 days. On the other hand, on short scales, 
(T < O.lTchar) it is the bottom topography that influences the most the model's 
solution. Both in the Black sea and in the square box the sensitivity to topog- 
raphy is approximately 40 times more important than the sensitivity to a. 

In addition to that, we should note that usually prescribed boundary condi- 
tions (impermeability and no-slip conditions have been used here as the initial 
guess for the cost function the minimization) seem not to be optimal for the 
model. As we can see in fig|T]and in figOl modifying a we can bring the model 
much closer to the high resolution model or to the observational data. Optimal 
numerical scheme allows to divide the cost function's value by 4 in the exper- 
iment with the Black sea model and by 3.5 in the experiment with the fourth 
order model in a square. But, the numerical scheme is modified in the assimila- 
tion process. As it has been shown in [18], optimal numerical scheme near the 
boundary may violate even impermeability condition indicating the necessity to 
change the domain's geometry. 

Taking into account an important influence of the numerical scheme that 
introduces boundary conditions into the model, it is reasonable to think about 
identification of the optimal scheme by data assimilation process instead of 
prescribing classical boundary conditions. 

Acknowledgments. The author thanks Gennady Korotaev from Marine 
Hydrophysical Institute, National Academy of Sciences of Ukraine for providing 
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21 



References 



[1] E. Lorenz, Deterministic non periodic flow., J. Atmos. Sci. 20 (1963) 130- 
141. 

[2] F.-X. Le Dimet, A general formalism of variational analysis, Tech. Rep. 
OK 73091, CIMMS report, Normann (1982). 

[3] F.-X. Le Dimet, O. Talagrand, Variational algorithm for analysis and as- 
similation of meteorological observations, theoretical aspects., Tellus 38A 
(1986) 97-110. 

[4] J.-L. Lions, Controle optimal de systemes gouvernes pas des equations aux 
derivees parielles., Dunod, 1968. 

[5] G. Marchuk, Formulation of theory of perturbations for complicated mod- 
els, Appl. Math. Optimization 2 (1975) 1-33. 

[6] M. Losch, C. Wunsch, Bottom topography as a control variable in an ocean 
model, J. Atmospheric and Oceanic Technology 20 (2003) 1685-1696. 

[7] A. W. Heemink, E. E. A. Mouthaana, M. R. T. Roesta, E. A. H. Volle- 
bregta, K. B. Robaczewskab, M. Verlaanb, Inverse 3d shallow water flow 
modelling of the continental shelf. Continental Shelf Research 22 (2002) 
465-484. 

[8] E. Kazantsev, 'Identifi cation of optimal topography by variational data assimilation!] 
J. Phys. Math. 1 (2009) 1-23. 



URL http : //www . dinaligroup . coin/volumel_2009_ JPM . htm 



[9] I. Shulman, Local data assimilation in specification of open boundary con- 
ditions, J. of Atmospheric and Oceanic Technology 14 (1997) 1409-1419. 

[10] I. Shulman, J. K. Lewis, A. F. Blumberg, B. N. Kim, Optimized boundary 
conditions and data assimilation with application to the m2 tide in the 
yellow sea, J. of Atmospheric and Oceanic Technology 15 (4) (1998) 1066- 
1071. 

[11] V. Taillandier, V. Echevin, L. Mortier, J.-L. Devenon, Controlling 
boundary conditions with a four-dimensional variational data-assimilation 
method in a non-stratified open coastal model. Ocean Dynamics 54 (2) 
(2004) 284-298. 

[12] P. G. J. ten Brummelhuis, A. W. Heemink, H. F. P. van den Boogaard, 
Identification of shallow sea models. International Journal for Numerical 
Methods in Fluids 17 (1993) 637-665. 

[13] F.-X. Le Dimet, M. Ouberdous, Retrieval of balanced fields: an optimal 
control method, Tellus A 45 (5) (1993) 449-461. 



22 



[14] Y. Leredde, J. M. Lellouche, J. L. Devenon, I. Dekeyser, On initial, bound- 
ary conditions and viscosity coefficient control for burgers' equation, Inter- 
national Journal for Numerical Methods in Fluids 28 (1) (1998) 113-128. 

[15] J. Lellouche, J. Devenon, I. Dekeyser, Boundary control on burgers equa- 
tion, a numerical approach, Comput. Math. Appl. 28 (1994) 33-44. 

[16] E. Kazantsev, Identification of an optimal boundary approximation by vari- 
ational data assimilation., J.Comp.Phys. 229 (2010) 256-275. 



[17] E. Kazantsev, Optimal boundary discretisation by variational data assimil ation"!] 
Int. J. for Numerical Methods in FluidsAccepted, doi: 10.1002/fld.2240. 
URL http : //hal . inr ia . f r/ inri a- 003888627^71 

[18] E. Kazantsev, Boundary conditions control for a shallow-water model.. Int. 
J. for Numerical Methods in FluidsAccepted. 

[19] X. Zou, I. Navon, F.-X. Le Dimet, An optimal nudging data assimilation 
scheme using parameter estimation, Q. J. of Roy. Met. Soc. 118 (1992) 
1163-1186. 

[20] V. Panchang, J. O'Brien, On the Determination of Hydraulic Model Param- 
eters Using the Strong Constraint Formulation Modeling Marine Systems, 
Vol. 1, CRC Press Inc, 1988, Ch. 1, pp. 5-18. 

[21] D. Chertok, R. Lardner, Variational data assimilation for a nonlinear hy- 
draulic model. Applied mathematical modelling 20 (9) (1996) 675-682. 

[22] A. Arakawa, V. Lamb, Computational Design of the Basic Dynamical Pro- 
cesses of the UCLA General Circulation Model, Vol. 17, Academic Press, 
1977, Ch. Methods in Computational Physics, pp. 174-267. 

[23] K. Ide, P. Courtier, M. Ghil, A. Lorenc, Unified notation for data assimila- 
tion: Operational, sequential and variational., J. of the Met. Soc. of Japan 
75 (IB) (1997) 181-189. 

[24] L. Hascoet, V. Pascual, Tapenade 2.1 user's guide[ Technical Report 0300, 
INRIA (2004). 



URL jhttp : //www . inria . f r/rrrt/rt-0300 . html 



[25] M.-H. Tber, L. Hascoet, A. Vidard, B. Dauvergne, 



Bui lding the tangent an d adj oint codes of the ocean general circulation model OPA with the automatic c 
Research Report 6372, INRIA (2007). 



URL https : //hal . inria. fr/inria-00192415 



[26] J. Gilbert, C. Lemarechal, Some numerical experiments with variable stor- 
age quasi-newton algorithms. Mathematical programming 45 (1989) 407- 
435. 

[27] C. Le Provost, J. Verron, Wind-driven ocean circulation transition to 
barotropic instability, Dyn.Atmos. Oceans 11 (1987) 175-201. 



23 



[28] E. Simmonct, M. Ghil, K. Idc, R. Temam, S. Wang, Low-frcqucncy vari- 
ability in shallow-water models of the wind-driven ocean circulation, J. 
Phys. Oc. 33 (2003) 729-752. 

[29] W. H. Munk, On the wind-driven ocean circulation, Journal of Meteorology 
7 (1950) 3-29. 

[30] G. Korotaev, T. Oguz, A. Nikiforov, C. Koblinsky, Seasonal, interannual, 
and mesoscale variability of the black sea upper layer circulation derived 
from altimeter data, JGR 108 (2003) 3122. 

[31] G. K. Korotaev, O. A. Saenko, C. J. Koblinsky, Satellite altimetry obser- 
vations of the black sea level, JGR 106 (2001) 917-934. 



24 



